PacBio HiFi 测序原理与应用¶

# ===== PacBio HiFi 数据处理流程 =====

# 1. 生成 HiFi reads（从原始 subreads 生成 CCS）
# 使用 PacBio 的 ccs 工具（SMRT Link 套件）
ccs input.subreads.bam output.hifi.bam \
    --min-rq 0.99 \       # 最低质量阈值 Q20（99%）
    --min-passes 3 \       # 至少 3 次完整 pass
    --num-threads 16

# 注意：Revio 平台直接输出 HiFi reads，无需此步骤

# 2. 查看 HiFi reads 统计
# 平均读长、总量、质量分布
python -c "
import pysam
bam = pysam.AlignmentFile('output.hifi.bam', 'rb', check_sq=False)
lengths = []
quals = []
for read in bam:
    lengths.append(read.query_length)
    quals.append(read.get_tag('rq'))  # 读质量
bam.close()
import numpy as np
print(f'HiFi reads 数: {len(lengths)}')
print(f'平均读长: {np.mean(lengths):.0f} bp')
print(f'中位读长: {np.median(lengths):.0f} bp')
print(f'N50 读长: {sorted(lengths, reverse=True)[len(lengths)//2]:.0f} bp')
print(f'平均质量: {np.mean(quals):.4f}')
print(f'总碱基数: {sum(lengths)/1e9:.1f} Gb')
"

# 3. 基因组比对
# 使用 minimap2（长读比对标准工具）
# -a: 输出 SAM 格式
# -x map-hifi: HiFi 专用预设参数
minimap2 -a -x map-hifi -t 16 \
    ref_genome.fa output.hifi.bam | \
    samtools sort -@ 8 -o aligned.sorted.bam
samtools index aligned.sorted.bam

# 4. 基因组从头组装（使用 hifiasm）
# hifiasm 是 HiFi 数据最佳组装工具
hifiasm -o assembly \
    -t 32 \                # 线程数
    output.hifi.fastq.gz   # HiFi reads

# 输出文件：
# assembly.bp.p_ctg.gfa → 主要 contigs（GFA 格式）
# assembly.bp.hap1.p_ctg.gfa → 单倍型1
# assembly.bp.hap2.p_ctg.gfa → 单倍型2

# GFA 转 FASTA
awk '/^S/{print ">"$2; print $3}' \
    assembly.bp.p_ctg.gfa > assembly.p_ctg.fa

# 5. 组装质量评估
# QUAST：基本统计
quast assembly.p_ctg.fa -r ref_genome.fa \
    -o quast_output/ -t 8

# BUSCO：基因完整性评估
busco -i assembly.p_ctg.fa -l mammalia_odb10 \
    -o busco_output -m genome -c 8